rm(list = ls())
options(width=110, device = 'quartz')

###################
## Assemble the Data 
####################
library(foreign)
library(plyr); library(dplyr)
library(coin)
library(rdrobust)
library(rms)
library(rdd)
library(ggplot2)
library(reshape2)
library(xtable)
library(Amelia)

POtoSouth11 <- function (pos, dta=st.info) {
  dta$South11[match(pos, dta$POAbrv)]
}
myRecode <- function(var, recode.ls, as.factor.result=TRUE, ...) {
  options(useFancyQuotes = FALSE)
  recode.ls <- sapply(recode.ls, function(x) paste(sQuote(x), collapse='='))
  recodes <- paste(recode.ls, collapse=';')
  print(recodes)
  revar <- Recode(var, recodes, as.factor.result=TRUE, ...)
  options(useFancyQuotes = TRUE)
  revar
}
rols <- function (..., data, subset=TRUE, cluster=1:nrow(data),
                  method="huber") {
  mod <- ols(..., data=data[subset, ], x=TRUE)
  robcov(mod, cluster=cluster[subset], method=method)
}
PasteEval <- function(..., print = FALSE) {
  txt <- paste(..., sep = "")
  if (print) { cat("\n", txt, "\n") }
  eval(parse(text = txt), envir = parent.frame())
}
FileName <- function (path=NULL, name=NULL, ext=NULL, replace=FALSE) {
  p <- paste(path, collapse = "")
  n <- paste(name, collapse = "")
  e <- paste(".", ext, sep = "")
  d <- format(Sys.Date(), "%y%m%d") 
  for (l in seq_along(letters)) {
    if (l > 1) old_file_name <- file_name
    file_name <- paste0(p, d, paste(n, letters[l], sep="-"), e)
    if (!any(grepl(file_name, list.files()))) {
      if (replace) file_name <- ifelse(l > 1, old_file_name, file_name)
      break
    }
  }
  cat("\nFile name:", file_name, "\n\n")
  return(file_name)
}
mysw <- function (...) {
  setwd(paste(..., sep = "/"))
}
mypdf <- function(name, ...) {
  date <- format(Sys.Date(), "%y%m%d")
  pdf(paste(paste(name, collapse=""), date, ".pdf", sep=""), ...)
}
tri <- function (x, h, c=0) pmax(0, 1 - abs((x - c) / h))


mysw(home.dir, 'Data')
load('party_control_data_161121.RData')
ls()
summary(data)


st.info <- read.csv("StateCodes.csv")


names(house <- read.dta('RD140410.dta'))
house <- plyr:::arrange(house, StPOAbrv, CDNumAtL, YearElec)
house <- group_by(house, StPOAbrv, CDNumAtL) %>%
    mutate(IncDWNOM1Lead1 = lead(IncDWNOM1, 1),
           IncDWNOM1Change1 = IncDWNOM1Lead1 - IncDWNOM1)

hnom <- read.dta('HL01112D21_PRES_BSSE.DTA')

################
#### Codebook
################
## demprop2: democratic proportion of the gubernatorial vote in year t
## Policy_lead_t1: gov’t Policy liberalism in year t + 1
## Policy_leads_delta2: change in gov’t liberalism between t and t+2

################
#### Setup RDD Variables
################

data.use <- subset(plyr::arrange(data, abb, year), !is.na(abb) & year >= 1936)

data.use <- data.use %>%
    group_by(
        abb
    ) %>% mutate(
        South11 = POtoSouth11(abb),
        DemGov = gov_party,
        DemPropGov0 = demprop2 - 1/2,
        DemMarginGov = DemPropGov0 * 2,
        DemWinGov = as.integer(DemPropGov0 > 0),
        DemSeatShareHouse = hs_dem_per_2pty / 100,
        DemSeatShareSenate = sen_dem_per_2pty / 100,
        DemControlHouse = hs_dem_control,
        DemControlSenate = sen_dem_control,
        GovLib = Policy,
        GovLibL1 = lag(GovLib, 1),
        GovLibL2 = lag(GovLib, 2),
        GovLibL3 = lag(GovLib, 3),
        GovLibL4 = lag(GovLib, 4),        
        GovLibL12 = (GovLibL1 + GovLibL2) / 2,
        GovLibP1 = lead(GovLib, 1),
        GovLibP2 = lead(GovLib, 2),
        GovLibP3 = lead(GovLib, 3),
        GovLibP4 = lead(GovLib, 4),
        GovLibP12 = (GovLibP1 + GovLibP2) / 2,
        GovLibD1 = GovLibP1 - GovLib, ## delta 1
        GovLibD2 = GovLibP2 - GovLib,  ## delta 2
        GovLibD3 = GovLibP3 - GovLib,  ## delta 3
        GovLibD4 = GovLibP4 - GovLib,  ## delta 4
        GovLibD12 = (GovLibD1 + GovLibD2)/2, ## ave of deltas 1-2
        GovLibDL1 = GovLib - GovLibL1 ## change, t relative to t-1
    ) %>% ungroup()
    
    


data.use <- group_by(data.use, abb) %>%
    mutate(GovLibL1 = lag(GovLib, 1),
           GovLibL2 = lag(GovLib, 2))

data.use <- mutate(data.use, GovLibL12 = (GovLibL1 + GovLibL2) / 2)

res.fe <- resid(lm(GovLib ~ abb + factor(year), data=data.use))
res.l1 <- resid(lm(GovLib ~ GovLibL1, data=data.use))
res.l2 <- resid(lm(GovLib ~ GovLibL2, data=data.use))
res.fe.l1 <- resid(lm(GovLib ~ abb + factor(year) + GovLibL1, data=data.use))
res.fe.l2 <- resid(lm(GovLib ~ abb + factor(year) + GovLibL2, data=data.use))

data.use$GovLib.fe0 <- NA
data.use$GovLib.fe0[as.integer(names(res.fe))] <- res.fe
data.use$GovLib.fe1 <- NA
data.use$GovLib.fe1[as.integer(names(res.fe.l1))] <- res.fe.l1
data.use$GovLib.fe2 <- NA
data.use$GovLib.fe2[as.integer(names(res.fe.l2))] <- res.fe.l2
data.use$GovLib.1 <- NA
data.use$GovLib.1[as.integer(names(res.l1))] <- res.l1
data.use$GovLib.2 <- NA
data.use$GovLib.2[as.integer(names(res.l2))] <- res.l2

data.use <- group_by(data.use, abb) %>%
    mutate(GovLib.1L1 = lag(GovLib.1, 1), ## lag of first-difference
           GovLib.1L2 = lag(GovLib.1, 2),
           GovLib.2L1 = lag(GovLib.2, 1), 
           GovLib.2L2 = lag(GovLib.2, 2),
           GovLib.fe0L1 = lag(GovLib.fe0, 1),
           GovLib.fe0L2 = lag(GovLib.fe0, 2),
           GovLib.fe1L1 = lag(GovLib.fe1, 1),
           GovLib.fe1L2 = lag(GovLib.fe1, 2),
           GovLib.fe2L1 = lag(GovLib.fe2, 1),
           GovLib.fe2L2 = lag(GovLib.fe2, 2),
           GovLib.1P1 = lead(GovLib.1, 1), ## lead of first-difference
           GovLib.1P2 = lead(GovLib.1, 2),
           GovLib.2P1 = lead(GovLib.2, 1), 
           GovLib.2P2 = lead(GovLib.2, 2),
           GovLib.fe0P1 = lead(GovLib.fe0, 1),
           GovLib.fe0P2 = lead(GovLib.fe0, 2),
           GovLib.fe1P1 = lead(GovLib.fe1, 1),
           GovLib.fe1P2 = lead(GovLib.fe1, 2),
           GovLib.fe2P1 = lead(GovLib.fe2, 1),
           GovLib.fe2P2 = lead(GovLib.fe2, 2),
           GovLibDL1 = GovLib - GovLibL1,
           GovLibP3 = lead(GovLib, 3),
           GovLibD3 = GovLibP3 - GovLib,  ## delta 3
           GovLibP4 = lead(GovLib, 4),
           GovLibD4 = GovLibP4 - GovLib)

################################################################################
#### ESTIMATES #################################################################
################################################################################


################################################################################
#### COMPARISONS ###############################################################
################################################################################

(house.cct <- with(house,
                   rdrobust(y=IncDWNOM1Change1, x=DifDPct, 
                            all=TRUE, kernel='tri')))

(nom.sds <- group_by(house, YearElec) %>%
     summarise(sd = sd(IncDWNOM1, na.rm=TRUE)))

house.cct$coef[1, ] /
    median(nom.sds$sd, na.rm = TRUE)

(pres <- subset(hnom, statenm == "USA" & cong >= 75))
pres$party2<-as.factor(pres$party)
pres.mod <- robcov(ols(dwnom1 ~party2, pres, x = TRUE),
                   cluster = pres$name)
                   

(pres.diff <- c(coef = coef(pres.mod)["party2=200"],
                se = sqrt(diag(vcov(pres.mod))["party2=200"])))

lib.sds <- group_by(data.use, year) %>%
    summarise(sd = sd(GovLib, na.rm=TRUE))
(med.lib.sd <- median(lib.sds$sd, na.rm = TRUE))
(med.nom.sd <- median(nom.sds$sd, na.rm = TRUE))

coef.df <- data.frame(Office = c(
                          "Governor\n(Policy)",
                          "State House\nMajority\n(Policy)",
                          "State Senate\nMajority\n(Policy)",
                          "U.S. House\nMember\n(Votes)",
                          "President\n(Votes)",
                          "State House\nMember\n(Votes)"),
                          ## "Governor\n",
                          ## "State House\nMajority\n",
                          ## "State Senate\nMajority\n",
                          ## "U.S. House\nMember\n",
                          ## "President\n",
                          ## "State House\nMember\n"),
                      Metric = c(rep("Policy", 3), rep("Roll Calls", 3)),
                      Estimate = c(.012, .029, .021, -house.cct$coef[3, ],
                          pres.diff[1], 1.186),
                      SE = c(0.006, 0.006, 0.004, house.cct$se[3, ],
                          pres.diff[2], 0.016),
                      sdxs = c(med.lib.sd, med.lib.sd, med.lib.sd,
                          med.nom.sd, med.nom.sd, 0.89))
(coef.df <- mutate(coef.df, StdEst = Estimate / sdxs, StdSE = SE / sdxs))


## mysw(home.dir, 'Graphics')
## mypdf("CoefCompare", width=7, height=4)
## (ggplot(coef.df, aes(x = reorder(Office, -StdEst),
##                      y = StdEst, ymin = StdEst - 1.96*StdSE,
##                      ymax = StdEst + 1.96*StdSE))
##  + geom_pointrange(size = .4)
##  ## + geom_hline(yintercept = 0, size = .3)
##  + labs(y = "Standardized Effect of Party Majority", x = element_blank())
##  + theme_bw()
##  + theme(panel.grid.minor = element_blank())
##  )
## dev.off()

mypdf("Figure8_CoefCompareLog", width=7, height=4)
(ggplot(coef.df, aes(x = reorder(Office, -StdEst), y = StdEst
                   ##, color = Metric
                   ## , ymin = StdEst - 1.96*StdSE, ymax = StdEst + 1.96*StdSE
                     ))
 + scale_y_log10(## breaks = c(.01, .1, 1)
                 )
 ## + geom_pointrange(size = .5)
 + geom_point(size = 4)
 + ylab(expression("Standaridized Party Effect (" * log[10] * " Scale)"))
 + labs(x = NULL)
 + theme_bw()
 + theme(panel.grid.minor = element_blank())
 + guides(color=FALSE)
 )
dev.off()

